Dynamics of a period-three pattern loaded Bose-Einstein condensate in an optical lattice 
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We discuss the dynamics of a Bose-Einstein condensate initially loaded into every third site of an optical 
lattice using a description based upon the discrete nonlinear Schrodinger equation. An analytic solution is 
developed for the case of a periodic initial condition and is compared with numerical simulations for more 
general initial configurations. We show that mean field effects in this system can cause macroscopic quantum 
self-trapping, a phenomenon already predicted for double well systems. In the presence of a uniform external 
potential, the atoms exhibit generalized Bloch oscillations which can be interpreted in terms of the interference 
of three different Bloch states. We also discuss how the momentum distribution of the system can be used as 
experimental signature of the macroscopic self trapping effect. 

PACS numbers: 03.75.F, 05.30.Jp 
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I. INTRODUCTION 

An optical lattice is a periodic potential formed by the ac- 
starkshift at the intersection region of two far detuned laser 
fields. A Bose-Einstein condensate loaded into such an opti- 
cal lattice is a useful system for studying and controlling the 
dynamics of ultra-cold atoms. In this system experimental 
studies in the meanfield regime have provided elegant demon- 
strations of band structure M |J, ^ and quantum chaos |Qj , and 
have considered the effects of nonlinearity on the condensate 
dynamics [|J |(| ffa. By suitably loading a condensate into deep 
optical lattices, number squeezing |Q| and the Mott-Insulator 
transition ||] have been observed. 

In this paper, we investigate a Bose-Einstein condensate 
loaded into every third site of an optical lattice, motivated by 
the recent experimental realization of this system by the NIST 
group [ |To| ] . In that experiment a combination of two indepen- 
dently controlled lattices was used to prepare the condensate 
into every third site of a single lattice. Briefly the procedure 
consists of loading a condensate into the ground band of lat- 
tice with periodicity 3a, so that the condensate is well local- 
ized in the potential minima of this lattice. A second lattice 
of periodicity a which is parallel to the first lattice, is then 
ramped up so that the superimposed light potentials form a 
super-lattice of period 3a. For the ideal case both lattice po- 
tentials are in-phase, and the addition of the second lattice will 
not shift the locations of the potential minima from those of 
the first lattice alone and the condensate will remain localized 
at these positions. Finally, by removing the first period-3a lat- 
tice on a time scale long compared to band excitations, but 
short compared to the characteristic time of transport within 
the lattice, the condensate will be left in every third site of the 
period-a lattice. 

A condensate loaded in this way is not an eigenstate of 
the final period-a lattice, and the condensate will continue 
to evolve in the final system. Outside the strongly correlated 
regime (where the depletion becomes large [|T"ij]) the Gross- 
Pitaevskii equation is expected to give a good description of 
the condensate dynamics, and has been applied to the lattice 
system 012, Ot|. In our system we assume the lattice is suf- 



ficiently deep that a tight binding description is applicable, 
and the Gross-Pitaevskii equation can be reduced to a discrete 
nonlinear Schrodinger equation involving only an onsite non- 
linear interaction and nearest neighbor tunneling. For the pe- 
riodic initial condition of equal wavefunction amplitude ev- 
ery third site and zero at all others, symmetry arguments can 
be used to reduce the wavefunction evolution to a two mode 
problem (analogous to a double well system with an energy 
offset) for which an analytic solution of the dynamics can be 
given. We show that for large ratios of the interatomic in- 
teraction strength to tunneling energy the condensate evolves 
with self-maintained population imbalance, whereby the con- 
densate population tends to remain localized in the initially 
occupied lattice sites. A similar phenomenon has been stud- 
ied in double wells systems [|l4[ 
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18[. We show 



that the momentum distribution of an interacting condensate 
changes in time in a manner which can be related to the spatial 
tunneling of condensate, and would be a suitable experimental 
observable. 

Under the influence of an external force a Bloch state will 
accelerate until it reaches the zone edge where it will Bragg 
scatter to the opposite side of the Brillouin zone. This peri- 
odic motion, known as Bloch oscillations, has been observed 
in systems of cold atoms JT^ ] and Bose-Einstein condensates 
in optical lattices [[}[ g], To illustrate how a linear potential 
affects the motion of a pattern loaded condensate we find ana- 
lytic solutions for the noninteracting case with periodic initial 
condition and show how the dynamics for this system can be 
interpreted in terms of the interference of three Bloch states of 
the lowest band Bloch oscillating in unison. We present nu- 
merical results for more general (non-periodic) initial condi- 
tions and consider characteristic properties of the momentum 
distribution. 

In Section H we start by reviewing the basic formalism of a 
Bose-Einstein condensate in an optical lattice in the tight bind- 
ing approximation. In Section |lj we study the dynamics for 
the case when only the optical potential is present. We discuss 
a periodic model and compare its predictions with numerical 
simulations for more general initial conditions. In section |y| 
we consider the dynamics of a Bose-Einstein condensate in 
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the lattice with the presence of constant external force. 



II. TIGHT BINDING DESCRIPTION OF A BEC IN AN 
OPTICAL LATTICE 



case the time scale for tunneling is determined by the hopping 
matrix element J. For this reason it is convenient to define a 
new scale of time t = Jt/h, and energies E n = e n /J and 
coupling constant A = NU / J. In terms of these new vari- 
ables the tight binding evolution equation (ph takes the form 



The dynamics of a condensate in the lattice can be modelled 
by the ID Gross-Pi taevskii equation 



ih 



V sm [ —x 



5$ h 2 <9 2 $ 

dt ~ ~2M'dx 2 

+y Ext (x)$ + iv^i|$| 2 $, 



(i) 



where V~Ext(x) is the external potential, Vq is the depth of the 
optical potential which is determined by the intensity of the 
laser beams, A is the wavelength of the lasers, a\ is the ID 
renormalized scattering length with dimension of an inverse 
length (which we obtain by requiring that the value of the 
chemical potential for the 3D system matches the ID chemical 
potential [pp||, M is the atomic mass and N the total number 
of atoms. 

In the tight-binding approximation, valid if the height of the 
lattice is higher than the chemical potential at each well, the 
ID condensate order parameter can be expanded in a Wannier 
basis keeping only the lowest band [^TJ : 



(2) 



= -(*n-l + *n+l) + (En + A|*„| 2 )*, 



(5) 



III. CASE OF NO EXTERNAL POTENTIAL 

A. The case of periodic initial conditions: reduction to a two 
mode system 

The stationary states and other aspects of the dynamics of 
Eq. (Q) have been discussed by others previously [^2, 23 , 24| |. 
We are motivated by recent experiments to treat the case of the 
particular nonstationary state described in the introduction, in 
which a condensate is initially loaded into every third site of 
an optical lattice. We treat first a model case in which no ex- 
ternal potential is present, (E n = 0), and in which the initial 
occupancies of each third site are the same, and in which the 
condensate initially has a uniform phase. At r = 0, the am- 
plitudes ^ n (r) are given by 



^3n(0) = y/p, 

*3»+l(0) - *3«+ 2 (0) = 0, 



(6) 
(7) 



where 4> n {x) = <p(x — na) is the condensate Wannier function 
centered on the n th lattice site, with J dx 4> n (x)4>* n+1 (x) = 
and J dx\4> n {x)\ 2 = 1, a = A/2 is the lattice constant, and 
^ n (t) is the n th amplitude. Replacing this ansatz in Eq. |lj 
the Gross-Pi taevskii equation reduces to the discrete nonlinear 
Schrodinger equation and the *„ satisfy: 



ih 



dt 



n+1) 



Nu\^ n \ 2 )^ n , (3) 



where U — {Aith? a\ j M) J efe|0 n | 4 is the strength of the 
on-site repulsion of two atoms occupying the lattice site 
n, e n = J VE X t(x)\4>n\ 2 dx describes the energy offset of 
each lattice site n, and J = / <j>* n+1 [-(h 2 /2M)(d 2 / dx 2 ) + 
Vq sin 2 (2ttx/ X)]<j) n dx is the hopping matrix element between 
adjacent sites. As shown in Ref. p2| , |5| , Eq. ([|) can be 
seen as an equation of motion ihd^ n /dt = dH /d^ 1 ^, de- 
rived from a Hamiltonian function H given by 

H = E - J (*n*«+l + K + l*n)+£n|*„| 2 + ^|*n| 4 , 

n 

(4) 

with i** and *„ treated as canonically conjugate variables. 
Both the Hamiltonian H and the norm J2 n \ ^n\ 2 = 1 are 
conserved quantities. 

The major interest of this paper is in the tunneling proper- 
ties of the condensate in the lattice, and in the noninteracting 



where Np is the initial number of atoms per occupied site. 
This initial condition is homogeneous in the sense that each 
occupied site has the same amplitude and phase along the 
length of the lattice. For an infinite lattice, or one with pe- 
riodic boundary conditions, the amplitudes for all initially oc- 
cupied sites (^>3n) evolve identically in time, and the ampli- 
tudes for the initially unoccupied sites satisfy ^3 n +i( T ) = 
^3n+2 (t) for all t and all n. This allows us to reduce the full 
set of equations (|^) to a set of two coupled equations 



= -2* 1+ A|* | 2 vI/ , 

OT 
.0*1 



8t 



-(tfi + tf ) + A|tfi| 2 #i, 



(8) 



(9) 



where *3 n = *o an d *3n+i = *3n+2 = *i for all n. The 
normalization condition is 

|* | 2 + 2|*i| 2 = p. (10) 

The Hamiltonian function, H, of this system is (from Eq. ^) 

H = -f-2(*S* 1 + * *J)-2|* 1 | 2 + ^|*o| 4 
p \ 2 



+A|*!| 4 

JAp 



(11) 
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By writing ^o^o-Jp and ^\=\j)\^ 
Eqs. (|)-(|) to the form 



d ( ipo 



dr \ ip 



I p/2, we can transform 



(12) 




where 7 = Ap is the ratio of on-site repulsion to tunneling 
energies, and the normalization condition ( |ic| ) is now \tpo\ 2 + 
\ipi\ 2 = 1. The factor of \/2 difference in the definition of ipo 
and ipi arises because ip\ represents the amplitude of the two 
initially unoccupied sites. With this factor incorporated, the 
matrix appearing in Eq. ( |l2| ) is explicitly Hermitian. We note 
this equation of motion is identical to that for a condensate in 
a double well trap in the two mode approximation [ ^5[ ^| . 

We note in passing that a similar reduction, to a system of 
[m/2] + 1 equations, exists for lattice systems that are loaded 
such that only every m-th site is initially occupied. 



B. Solution of the equations of motion 

It is convenient to write the lattice site amplitudes appearing 
in Eqs. (|) and @ as ^ =fe ie ° y /p and * x = ge i9l jp, where 
/, g, and 9 are real. By introducing the phase difference <f> = 
6*0-6*1, Eqs. (§), (ph and dHl) can be recast as 



/ = 

9 = 

1 
2 



2g sin0, 
-fsintj), 

-Afgcos4>+l(f i + 2g i )-2g 2 



(13) 
(14) 

(15) 



The analytic solutions of Eqs. (|13j)-(|15j), found using a pro- 
cedure similar to that presented by Raghavan et al. Jl4|], 
can be expressed in terms of Weierstrassian elliptic functions 
p(r; c/2, 93) [p7|- The analytic solutions are 



f(r) 



9{r) 



24 



12p(T;<7 2 ,33) + (9 + 2 7 + 7 2 ) 



(16) 



(17) 



where the parameters g 2 and g 3 are given by 

.92 = (81 - 14 7 2 + 4 7 3 + 7 4 )/12, 

(729 + 243 7 2 - 46 7 3 - 15 7 4 + 7 6 ) /216. (19) 



(18) 
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The solutions /(r) and g(r) are oscillatory functions whose 
amplitudes and common period, T( 7 ), are determined by the 
parameter 7 (see Figs. 1 and 2). It is useful to qualitatively 
divide this behavior into two regimes, separated by 7 = 2. 
Analysis of Eqs. (|3j) - (|5j) shows that /(to) = <?(to) for 
some value of To when 7 < 2, and for 7 > 2, /(r) > <?(t) for 
all t. 
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FIG. 1: Oscillation period (in units of ft/ J) as a function of the in- 
teraction strength 7. 



In this case the role of interactions is relatively small, 
and the behavior can be approximately understood by taking 
7 = 0, in which case the matrix of Eq. ( p"2| ) is constant in 
time. The equations of motion in this case are equivalent to 
those of a two-state Rabi problem where the two levels 
are coupled by a Rabi frequency of strength y/2, which is de- 
tuned from resonance by —1. This system will undergo Rabi 
oscillations whereby atoms periodically tunnel from the ini- 
tially occupied site into the two neighboring sites. Because 
the coupling is detuned from resonance the transfer of pop- 
ulations between wells is incomplete, with \ipi\ 2 attaining a 
maximum value of 8/9. The Rabi model predicts that the cy- 
cling frequency between the levels is equal to the difference 
between the eigenvalues of the matrix of matrix of Eq. (|l2|), 
which gives the period of oscillation as TR a bi = 2ir/3 in units 
of hj J (see Fig. [j}. 



2. Interaction dominated regime 

The effect of interactions on the mean-field dynamics is to 
cause the energies of the initially occupied sites to shift rel- 
ative to those of the unoccupied sites. As 7 increases and 
this energy shift increases relative to the strength of coupling 
between sites, the tunneling between sites occurs at a higher 
frequency, but with reduced amplitude. The population of the 
initially occupied sites becomes self trapped by the purely re- 
pulsive pair interaction, which in the context of a double well 
system has been called "macroscopic quantum self trapping" 
P6|]. This is demonstrated quantitatively in Fig. ^ where we 
plot the minimum value of f 2 occurring during the oscilla- 
tion as a function of 7. In contrast to the tunneling domi- 
nated regime, where tunneling periodically populates all sites 
equally, the condensate tends to be localized on the initially 
occupied sites in the interaction-dominated regime. 



1. The tunneling dominated regime 



C. Momentum space dynamics 



For 7 < 2, we find that the oscillation period is essentially 
constant (see Fig. 1). 



Typically the spacing between individual wells in an opti- 
cal lattice is too small to resolve the localized density distri- 
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FIG. 2: Minimum value of f 2 during an oscillation period as a func- 
tion of 7. As 7 increases the population imbalance between wells 
increases (see text). 



butions of atoms in neighboring sites using standard imaging 
techniques. 

The momentum distribution is a more convenient observ- 
able which approximately corresponds to the expanded spa- 
tial distribution of the released condensate. Here we calculate 
the momentum dynamics of the condensate loaded into ev- 
ery third site of an optical lattice, and show how this relates 
to evolution of the spatial amplitudes given in Eqs. (16) and 

(ITT 



J). 

In the tight binding approximation the condensate order pa- 
rameter $(2;, r) is expressed as a sum over the lattice sites (g). 
Because of the periodicity of the system, the momentum space 
wavefunction, which we denote as $(fc, r), is expressible as 
the Fourier series 



$(fc,r) = ^* m (r)e^ m X (fc), 



where 



e ikx (j)o{x)dx. 



(20) 



(21) 



To compute the momentum distribution, we invoke the 
identity J2n=~o^ ^ kna = N s 6 k ^ m/a , where N s is the num- 
ber of lattice sites, and m an integer. Since there are only two 
independent amplitudes in the set {'I'm}, we find that 



$(fc,r) = \/3jpc rn (T)x m Sk, qm /3, 



(22) 



(r) = (* + 4-ie-' 9 " 1 / 3 + m 2 e- l2qm / 3 ^j 



where q is the reciprocal lattice vector q = 2n/a. The mo- 
mentum distribution of the system has very sharp peaks of rel- 



ative amplitude 



at momentum k = qm/3, arising from 



the 3-lattice site spatial periodicity of the condensate wave- 
function. In addition, \ m describes a slowly varying envelope 
determined by the localization of the Wannier states at each 
lattice site. 



Using the analytic solutions for / and g, and Eq. (|l5|), we 
obtain only two independent Fourier amplitudes 



|c 3 n(r)| 2 

|c 3 n+l(T)| 2 



|c 3 n+2(r)| 5 



(3/ 2 + l)(/ 2 -l)J, (23) 

(3/ 2 + l)(/ 2 -l)) (24) 
(25) 



In the reduced zone scheme, where we only consider momenta 
in the range k <E [— q/2, q/2], the momentum wavefunction 
then consists of three peaks corresponding to Bloch states of 
quasimomenta 0,±q/3. The identical behavior of |c3„+i| 2 
and |c3, i+ 2 1 2 means that the ±g/3 peaks always have the same 
intensity. If interactions between the atoms are ignored (i.e. 
7 = 0), the momentum components are constant in time (see 
Eqs. ([23|) and (^)), even though tunneling occurs between the 
lattice sites. However, when interactions are considered, the 
momentum intensities explicitly depend on the occupations 
of each site and will vary in time when tunneling occurs. The 



magnitude of the time variation of the 



is proportional to 



7, but will reduce for sufficiently large values of 7, where the 
self-trapping effect causes the tunneling between lattice sites 
to stop (i.e. f 2 w 1 at all times). In Fig. || we show the 
maximum contrast between the intensity of the Fourier peaks, 
defined as A max = (|ci| 2 — |co| 2 ) max , where the value is 
maximized by evaluating the |c„| 2 at the time when / takes 
its minimum value (we note that /(r) is given by Eq. (|T6|)). 
We note that the maximum contrast occurs for 7 w 3.9; in 
this case, the zero quasi-momentum component cq(t) van- 
ishes once during each period of oscillation T(j). For values 
of 7 greater than 3.9 the contrast between the intensities starts 
to decrease due to the reduction in tunneling caused by the 
nonlinearity- induced self-trapping. 
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FIG. 3: Maximum contrast of the Fourier components as a func- 
tion of 7. The maximum contrast is defined as A max = (|ci| 2 — 
|co| 2 )max with the maximum value occurring when f 2 is at its min- 
ima. 



D. Application to an inhomogeneous condensate 

Here we wish to consider the dynamics for an inhomoge- 
neous condensate, applicable to a condensate initially pre- 
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pared in a harmonic trap. For the pattern loaded condensate, 
we use inhomogeneous to refer to the overall spatial envelope 
of the period three initial condition. The previous homoge- 
neous theory we have presented is expected to accurately de- 
scribe inhomogeneous cases when the initial pattern of pop- 
ulation in every third site extends over many lattice sites i.e. 
]V S > 1 so that mean-field energy associated with each triplet 
of sites U(n) = 4 Xa=i l^3n+j| 4 varies slowly across the 
system. Taking as a particular example we choose a gaussian 
envelope to the periodic arrangement of atoms into every third 
site, so that the initial state is 



*3n+l(0) = *3n+ 2 (0) - 0, 



(26) 



*3n(0) = 3 



rJV? 



1/4 



exp(-(9n/7V s ) 2 ). (27) 



For the simulations shown here, the parameters used were 
N = 10 5 , N a = 76, U = 2.11 x KT 5 £ R and J = 0.075-Br, 
where E R = h 2 /8Ma 2 is the lattice recoil energy; these cor- 
respond to a condensate of 10 5 atoms of 87 Rb produced in a 
magnetic trap with axial and radial frequencies of 9 and 12 Hz 
respectively, and loaded into a lattice of counter-propagating 
light 785nm of depth 4.5-Er, with E R = 2.2kHz. These pa- 
rameters are typical of the experimental regime, but also lie in 
a range in which the homogeneous model is expected to give 
a fair description. The total number of occupied wells and 
the strength of the on site interatomic interaction were cal- 
culated by preserving the value of chemical potential of the 
system upon reduction to one spatial dimension and by as- 
suming that each of the localized orbitals in the tight binding 
description are Gaussian. The hoppingrate J was estimated 
by using Mathieu functions. In Fig. |] we show the evolu- 
tion of the population of the central wells (normalized to one) 
compared with the homogeneous model with 7 taken to be the 

local mean-field energy 7 e g = Ay-| (9/N s ). 

Figure [|, shows the results of numerical integration of the 
equations of motion and the approximate analytical solution 
given by the quasi-homogeneous model described above. We 
see that the numerical and analytical results agree well at short 
times, but differ more as time progresses due to the different 
mean-field seen by different wells. 

To understand the disagreement as time evolves, we show in 
Fig. ||the numerical Fourier spectrum for the inhomogeneous 
case evaluated at several different times. The variation of the 
intensities of the peaks, which, as shown below, is related to 
the spatial tunneling between lattice sites in the presence of 
the mean field, can be seen in the plot. Initially all occupied 
sites are in phase and the three distinctive momentum peaks 
have a narrow width determined by the intrinsic momentum 
uncertainty of the condensate envelope. That is the reason why 
the homogeneous model fits very well. However, as time pro- 
gresses the meanfield variation across the lattice causes the 
tunneling rates to vary with position and leads to the momen- 
tum peaks broadening. The homogeneous description then, 
starts not to be very accurate. 

We note that momentum space signature for spatially tun- 
neling in the interacting system is still present in the inho- 
mogeneous case. This is shown in Fig. ^, where we plot the 
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FIG. 4: Comparison between the population evolution of the central 
three wells for the inhomogeneous condensate and the homogeneous 
model, inhomogeneous condensate: stars for the initially populated 
well and boxes for the initial empty wells. Homogeneous model: 
dashed line represents the initially populated wells, and the solid line 
represented the initially unpopulated wells. We used jgg as the local 
mean field energy. The parameters used for the simulation were J = 
0.075.Z5R, and y e g = 2.64. The time is in computational units: t = 
tsih/En where En is the recoil energy and tsi is time in SI units. 
(Wh/E R « 0.7ms). 



t =0 



t = 12 




q/3 2q/3 
momentum 



q73 2q73~ 
momentum 



FIG. 5: Momentum distribution of the inhomogeneous condensate 
evaluated at various times for the same parameters as used in FigM 



evolution of the quantities |c 3 „(t)| 2 , |c 3ii+ i(t)| 2 , |c 3n+2 (T)| 2 
(calculated from the numerical simulation by partitioning the 
numerical Fourier spectrum in three equal non overlapping 
sections, each centered around the respective peak and adding 
the square of the norm of the Fourier components within each 
section) vs. the ones calculated with the homogeneous model, 
but using an averaged value 7 avc = |^n| 4 instead of 

7 = Ap. It can be observed that the predictions of the simple 
model are in very good agreement with the numerical results 
when the three peaks of the spectrum are well defined. For 
longer times, the width of the Fourier peaks increases, until a 
point when they split. At this point the quantities |c 3 „(t)| 2 , 
l c 3n+i( r )| 2 > l c 3n+2(' r )| 2 are not meaningful anymore. Be- 
cause the parameters used for the numerical calculations were 
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chosen to be experimentally achievable, and as shown in the 
plots the model predictions are fair at least for one period of 
oscillation, we conclude that Fourier distribution can be used 
as a signature of the mean-field quantum tunneling inhibition. 
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FIG. 6: Evolution of momentum peak populations. Upper curves: 
population of the q = ±27r/3 momentum states. Lower curves: 
population of the q — momentum state. Inhomogeneous con- 
densate (dotted), homogeneous result (solid line), where the com- 
parison is made by replacing 7 by an average mean field energy 
. = A^„ \^ n \ = 1.85. Parameters are the same as in Fig. 



IV. DYNAMICS WITH A CONSTANT EXTERNAL FORCE 

A. Homogeneous three state model 

In this section we consider the dynamics of a periodically 
loaded condensate in the presence of a linear external poten- 
tial, corresponding to a uniform force parallel to the lattice. 
In what follows we assume that the force is sufficiently weak 
that band excitations due to Landau-Zener tunneling is neg- 
ligible, so that a tight binding picture of the lowest band is 
sufficient to describe the dynamics. In this case the evolution 
equation differs from what we considered in the previous sec- 
tion by the term E„ in Eq. (Q) taking the form E n = ri£, 
where £ is the potential difference between lattice sites (in 
units of the hopping matrix element J). Taking the initial 
conditions (||)-(7), and transforming the Wannier amplitudes 
as * 3 „+,(t) = V 3n+J (T)e- l3 < T (j = 0, 1, 2), we obtain the 
evolution equations 



1 dr 


= -(* 3 „_ ie l3 « T + * 3 „ +1 ) 






+A|^ 3n | 2 ^ 3n , 


(28) 


% dr 


= -(*3n+2 + *3n)+£*3fM 


-1 




+A|*3n+l| 2 *3n+l, 


(29) 


.3§3n+2 

1 dr 


= -(f 3n+ i + f 3 n+3e" i3Cr ) 


+ 2^3n+l 




+A|* 3 „ + i| 2 * 3n+1 , 


(30) 



Assuming periodic boundary conditions, the periodicity of the 
initial conditions and equations of evolution allow consider- 
able simplification from the full set of 3T coupled equations. 
In particular these assumptions mean that every third Wannier 
amplitude evolves identically (i.e. & n = ^ n +3) and so the 
evolution of the system can hence be reduced to the three in- 
dependent equations 



l ~dr- 

.8^2 



-(* 2 e i3?T + *i)+A|* | 2 *o, (31) 

-(* 2 + *o)+£*i+A|*i| 2 *i, (32) 

(* 1 +* e- !3 « T ) + 2C* 2 

-A|* 2 | 2 * 2 , (33) 



where the new amplitudes map onto the original set according 
to * «-> {* 3 „}, *! «-> {*3n+i}, and § 2 <-> {tf 3n+2 }, and 
obey the normalization condition 



(34) 



The equations of motion (|31[)-(|33[) are more difficult to treat 
analytically than the case considered in the last section due to 
the presence of a linear potential. In this paper we derive an 
analytic solution for a noninteracting condensate (i.e. A=0), 
which provides valuable insight into the complicated tunnel- 
ing dynamics the system exhibits in the absence of nonlinear- 
ity, yet should furnish a good description for dilute conden- 
sates satisfying 7 < 1. For the nonlinear regime we present 
numerical results to illustrate typical behavior. 



B. Analytic solution for linear dynamics 

Defining the vector x(i) = (*o(r), *i (r), ^(r)), and us- 
ing the transformation y(r) = P(t)x(t), where P(r) is the 
unitary matrix 



1 1 1 

P( T ) = I e~ ifT e - '^-^ e-'C^+^f) 

-i2£r e - I (25r+^) e _i(2£r-4f) 



(35) 



the linear version of Eqs. (|3jJ)-([33|) can be decoupled, directly 
yielding the solutions 



3VT 

where we have defined the phase terms 



^ e -iAo(r) +e i(a 3 a_A 1 (r)) 



A n (r) 
forn = 0, 1,2. 



2mr 2nir 
sm(£r — ) +sin(— ) 



(36) 



(37) 
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These solutions for the spatial amplitudes ^(r) can be 
most easily understood by considering a Bloch state decom- 
position of the condensate wavefunction. The nature of our 
system allows us to construct an analytic form for the initial 
wavefunction. Because the system has a three lattice site pe- 
riod and is assumed to be in the lowest band, the wavefunction 
can beO expressed as a superposition of three Bloch waves (of 
the lowest band) which are symmetrically spaced in quasimo- 
mentum. Assuming the condensate initially has a total crystal 
momentum of zero, at this time the wavefunction must be of 
the form 

$(z, 0) = a <p (x) + a+ip q/3 (x) + a_ip_ q/3 (x) (38) 

where ipk( x ) is a Bloch state with quasimomentum k, and the 
a are complex constants determined by the lattice depth, with 
|a+| = \a-\. 

The action of an external force on a Bloch state causes it to 
linearly change its quasimomentum in time according to 



T + fc(0). 



(39) 



The periodicity of the Bloch dispersion relation in k, and 
hence of the group velocity of the Bloch wave, gives rise to 
the well-known phenomenon of Bloch oscillations [^9|. For 
the case we are considering here, the system consists of three 
Bloch states whose quasimomenta will translate in unison un- 
der the action of the external force. During this evolution each 
state accumulates phase at a rate determined by the instanta- 
neous Bloch energy, i.e. 



A n (r) 



E(k n (s))ds, 



(40) 



where fc„(r) = — £r/a + k n (Q) is the quasimomentum of 
Bloch state n at time t. In the tight binding approximation the 
dispersion relation for the Bloch states has the analytic form 



E(k) = -2cos(fca), 



(41) 



for which A„(r) can be evaluated,and yields the results given 
in Eqs. (p7|)). The wavefunction evolution in the Bloch basis 



+a-ip_t t _ i {x)e- lA ^ T \ (42) 

From this solution we can obtain solutions for the evolution 
of the spatial amplitudes Eqs. ([37]). To do this we expand 
the Bloch states in terms of Wannier functions according to 



- lkna 4>n( x ) and make use of Eq. (g). Note: 
we take all a = | as determined by the initial conditions, 



<Pnk(x) = £„ 



Eqs.(|) and (@). 



1. Bloch Oscillations 



The evolution of the spatial amplitudes, and in particular the 
population in each well is then determined by the interference 



of the Bloch phases A„. These functions are all periodic in 
time with period tb = 27r/£ (in units of r = tJ/K). This is 
the normal Bloch oscillation period, and gives the time scale 
over which the quasimomenta of the Bloch states develop by 
exactly a lattice vector. 



2. Small £ solution - Non classical transport 

To understand the dynamics, we first start by considering 
the case when £ is small. For this case, the population in the 
wells is given by 



l*o(r)| 2 
l*i(r)| 2 
l* 2 (r)| 2 



1 



5 + 4cos(3t) +£ 2 1i{t) 



(43) 



^-((2 + 30 [1 - cos(3r) - ^h(r)])44) 
((2 - 3£) [1 - cos(3r) - |V)])45) 



where 



h(r) 



T 

T 



3r + 3tcos(3t) — 4sin(3r)) . 



(46) 



The above solution shows that when the force is applied, the 
degeneracy in the population of the wells represented by ^ 
and ^2 is lifted. For £ > 0, atoms in Wo start to tunnel to 
^fi more rapidly than to fy 2 - This should be compared with 
the results in the absence of the force, where Wi and ^> 2 be- 
have identically. Thus, in this weak limit, the effect of a linear 
potential to enhance the tunneling from the initial populated 
3n wells to their in + 1 neighbors ones, making the system 
closer to resonance, in the sense of Eq. ( 12), where the reso- 



nance condition results in the initially populated wells becom- 
ing empty at some later time. It is interesting to note that the 
system exhibits "nonclassical" dynamics whereby the atoms 
start to tunnel in the direction opposite the direction of the 
force (this statement applies even when the external field is 
not weak). 



3. Resonances 

In Fig. ^ we show the temporal evolution of the spatial 
amplitudes ^ for a range of values of £. For certain choices of 
£ the initially occupied \Po amplitude periodically disappears 
- we refer to these as resonances. By requiring = in 
Eq. (^6j) we obtain the following conditions on the phases for 
these resonances 



Ai(r) - A (r) = ± (~ + 2 7 



A 2 (t) - Ai(t) = ± 



2-Km 



(47) 
(48) 



where n and m are integers and the same sign choice on 
the right hand side must be made for both equations. When 
these conditions are satisfied the population is not shared be- 
tween the adjacent sites, but preferentially tunnels to one of 
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FIG. 7: Evolution of the normalized population for different values 
of £. One Bloch period is shown in the plots except £ = where the 
period is infinite. Solid line: l^nl 2 , dotted line: I'^n+i j 2 , dashed 
line: |^3 n +2| 2 . The "nonclassical" motion can be seen where the 
3n+l-well populations increase more rapidly that the populations of 
the 3n+2-wells. It can also be seen in the plots that £ = £ max /2 and 
£ = Cmax are resonant values. 



which is shown in Fig. ^. For large \n\ or |m| the resonant 
values of £ are close together, and become spaced further apart 
as the magnitude of n decreases. 

In the absence of an applied force (i.e. £ = 0) the phase 
terms are time-independent with Ao = — 2, Ai = —1 and 
A 2 = —1 so that the resonance conditions can never be 
achieved. When the external force is applied the phases os- 
cillate at a rate that increases with £ and an amplitude that 
decreases with £. The largest value of £ for which a reso- 
nance can be found is £, r °|x = 6Vo/27r, since for values of £ 
greater than this the amplitudes of the phase oscillations A< 
are so small that they cannot satisfy the condition for popu- 
lation resonance. In the regime £ > £,™l K , | ^0 1 exhibits only 
one (nonzero) minima per Bloch period and the dynamics of 
the system is dominated by the Bloch oscillations. 



4. Large Force Limit 



For £ > £, r °a X the system exhibits a population imbalance 
as the Bloch oscillation suppresses the ability of the system to 
tunnel between sites. In the limit £ 3> tne population in 
the wells is described by rapid, small amplitude oscillations 
around its initial value, 



l*o(r)| 2 
l^iWI 2 



3 



1 8 • 2 & 
'--Bin \- 



12 1 

w s e 8m ' 



(50) 
(51) 
(52) 



j- res 

6 V3 



The case of £ = 3£^x depicted in Fig. 8 shows an approach 
to this behavior. 



C. Numerical results for the non-linear case 



n 



FIG. 8: The spectrum of values of external force, ordered in decreas- 
ing magnitude, for which a population resonance occurs i.e. the val- 
ues of £ for which >i?o periodically disappears 



the neighbors. In particular the +sign choice in Eqs. ( |47| ) and 
d48|), gives the phase condition for complete tunneling into 
^1. Similarly the —sign case gives the condition for complete 
tunneling into ^2- Solving for the values of £ for which Eqs. 
( p^ ) and ( p^ ) hold, gives the spectrum 



C res 



± 



6^3 



(1 + 3n) 



7T (l + 3n) 2 + 3(n + 2m+l) 2 



(49) 



Because of the difficulty of solving analytically the equa- 
tions of motion the nonlinear equations were solved numer- 
ically. Although the dynamics in this situation is consider- 
ably more complicated than in the linear case, the resonance 
picture still gives us useful guidance concerning the expected 
behavior. In general, there exist critical values of £ and A 
above which no resonances occur. Furthermore, nonlinearity 
tends to destroy the periodicity of the Bloch oscillations that 
is present in the noninteracting case. For example, in Fig. ^|, 
we see that for £ = 2£™| x , the introduction of interactions 
brings the system into resonance, but that for large interaction 
strength, nonresonant behavior is restored. For£ = 2££^ x /7, 
on the other hand, the introduction of interactions eventually 
draws the system out of resonance. 

Concerning the momentum distribution, similarly to the un- 
tilted case, interatomic interactions induce time variation of 
the momentum intensities. The contrast between momentum 
components vanishes at A = 0. As A increases, the contrast 
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FIG. 9: Effects of interactions on generalized Bloch oscillations for 
the pattern loaded system. Evolution of \SS? n 2 for various interaction 
strengths. Upper plot: £ = 2£S X . Lower plot: £ = 2£^ x /7. 



increases to a maximum value, then eventually decreases to- 
wards zero when macroscopic imbalance is achieved, analo- 
gous to what is seen in Fig. |3| for the untilted case. Because 
the external field together with the nonlinearity breaks down 
the periodicity of time evolution, the dynamics of the momen- 
tum distribution is quite complex. 

To compare the predictions of the homogeneous model in 
the presence of a linear external potential to a more realis- 
tic case, we again solved numerically the discrete nonlinear 
Schrodinger equation for a condensate loaded every three lat- 
tice sites but, instead of being homogeneous, initially with a 
Gaussian profile. Very good agreement between the model 
and the numerical results was found for short times and mod- 
est mean-field energies, if, as in the untilted case, we use an 
effective mean field energy. In Figs. [k] and 1 1 we present a 
comparison between the evolution of the normalized popula- 
tion at the central wells found numerically and the prediction 
of the model for the parameters: NU = 4.8Sr, N s = 290, 
J = 0.075-Er and £ = 2 with computational time of 50 is 
about 3.5 ms. For longer times, the model predictions start to 
disagree with numerical simulation due to the spread of the 
three Fourier peaks. We observe that the dynamics are much 
more sensitive to inhomogeneous effects in the presence of an 
external force, and for the more inhomogeneous initial state 
used in Figs. 4-6 (which extends over N s = 76 sites), the in- 
homogeneous result much more rapidly departs from the ho- 
mogeneous prediction than the example presented here. 




FIG. 10: Comparison between the evolution of a inhomogeneous 
condensate with the homogeneous result. Inhomogeneous conden- 
sate: |* | 2 (boxes), |*i| 2 (stars), and |* 2 | 2 (diamonds). Homo- 
geneous case: l^nl 2 (dash-dot line), l^n-i-il 2 (solid line), and 
|^3n+2| 2 (dotted line), where we have taken 7 as the local mean 

2 and 



field energy. The parameters used were J 
Teff = 1-59 
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FIG. 11: Evolution of the momentum peak populations, q = 
(dashed line, squares), q — —2n/3 (dash-dot line, stars), q — 2n/3 
(solid line, diamonds). Inhomogeneous condensate (dotted curves), 
homogeneous model (lines) using the same parameters as those in 
Fig. llQ, but replacing 7 by an average value 7 avc = 1.11 for the 
homogeneous model. 



V. CONCLUSION 

We have studied the dynamics of a condensate loaded in an 
optical lattice in such a way that initially only every third well 
in the lattice was occupied. We considered two cases, one 
with only the optical potential, and one with a superimposed 
constant external field. It was found that in both cases inter- 
atomic interactions can cause macroscopic quantum self trap- 
ping, which is a self maintained population imbalance across 
the wells. The analysis was based on a tight binding model 
that assumes that at t = all the 3n wells were identically 
populated. Using this model, we studied the different regimes 
where the self trapping phenomenon occurs and estimated the 
critical parameters for each case. We have developed an ana- 
lytic solution for the dynamics of a periodic initial condition in 
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the absence of an external potential. For this case the symme- 
try of the system reduces the discrete nonlinear Schrodinger 
equation to a two-mode problem with an elliptic function so- 
lution. We have verified the usefulness of the analytic solution 
by comparing it to numerical solutions of the discrete nonlin- 
ear Schrodinger equation for more general initial conditions. 

We have demonstrated that meanfield effects cause the mo- 
mentum distribution to vary with time, and have shown how 
this variation relates to the spatial tunneling in the system. 
This result suggests that the temporal variation in the time-of- 
flight density distribution (which approximately corresponds 
to the in-situ momentum distribution) would be a useful ex- 
perimental signature of the spatial tunneling. 

In the presence of an external potential and neglecting in- 
teractions, it was found that the spatial population amplitudes 



evolve periodically in time, with the oscillation period deter- 
mined by the force parameter £. We have shown that these 
dynamics can be understood as the interference of three quasi- 
momentum states Bloch oscillating in unison. When interac- 
tions are taken into account, the periodicity of the system is 
destroyed and the atoms exhibit more complicated dynamics. 
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